1  Prognosemodel Retentie na 1 jaar

ITD | B Communication and Multimedia Design (CMD) - voltijd - versie 1.0

Auteur

Theo Bakker, lector Learning Technology & Analytics, De HHs

Publicatiedatum

20 januari 2025

1.1 Methode, data en analyse

1.1.1 Toelichting op de methode

Voor de ontwikkeling van prognosemodellen gebruiken we de aanpak van Tidymodels. Tidymodels is een framework voor het bouwen van een prognosemodel. Hiermee verzekeren we ons van een systematische, herhaalbare en schaalbare aanpak.

1.1.2 Toelichting op de data

De basis voor deze analyse is studiedata van De Haagse Hogeschool (De HHs), verrijkt door het lectoraat LTA. De data bevat informatie over de inschrijvingen van studenten in het eerste jaar van de opleiding:

  1. Demografische kenmerken: geslacht, leeftijd, reistijd en SES totaalscore.
  2. Vooropleidingskenmerken: toelaatgevende vooropleiding, studiekeuzeprofiel en gemiddeld eindcijfer in de vooropleiding.
  3. Aanmeldingskenmerken: aansluiting (direct na diploma, tussenjaar, switch), dag van aanmelding, aantal parallelle studies aan De HHs en collegejaar.

Deze variabelen zijn gekozen omdat we uit eerder onderzoek weten dat ze voorspellende waarde kunnen hebben voor studiesucces (Bakker, 2022) of omdat ze behoren tot sensitieve kenmerken die in fairness analyse gebruikt worden.

Toon code
# Read the data dictionary
dfData_dictionary <- Get_Data_Dictionary()

# Show the data dictionary
Get_tblData_Dictionary(dfData_dictionary)
Tabel 1.1: Variabelen en mogelijke waarden

Variabele

Toelichting

Mogelijke waarden

Aanmelding

De dag van de aanmelding voor de studie gerekend vanaf 1 september

-366 - 366

Aansluiting

De manier waarop een student instroomt in de opleiding

2e studie, Direct, Na CD, Onbekend, Overig, Switch Extern, Switch Intern, Tussenjaar

APCG

Of de student ten tijde van het behalen van de toelaatgevende vooropleiding in een armoedeprobleemcumulatiegebied woonde

Ja, Nee, Onbekend

Cijfer_CE_Engels

Het gemiddelde cijfer Engels van het centaal examen van de middelbare school

0-10

Cijfer_CE_Natuurkunde

Het gemiddelde cijfer voor Natuurkunde van het centaal examen van de middelbare school

0-10

Cijfer_CE_Nederlands

Het gemiddelde cijfer Nederlands van het centaal examen van de middelbare school

0-10

Cijfer_CE_VO

Het gemiddelde cijfer voor het centaal examen van de middelbare school

0-10

Cijfer_CE_Wiskunde

Het gemiddelde cijfer voor Wiskunde van het centaal examen van de middelbare school

0-10

Cijfer_SE_VO

Het gemiddelde cijfer voor het schoolexamen van de middelbare school

0-10

Collegejaar

Het collegejaar van de inschrijving

2012-2022

Dubbele_studie

Gegeven of de student meer dan 1 studie volgt

Ja, Nee

Geslacht

Het geslacht van de student

M, V

ID

Uniek nummer per student

Leeftijd

De leeftijd van de student op 1 oktober

16-100

Ranking

De positie die de student heeft in de selectie voor de opleiding (alleen bij de B Huidtherapie)

0-350

Reistijd

De reistijd van de student naar de vestiging van de opleiding op een maandag om 9 uur met het openbaar vervoer vanaf de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding in minuten

0-120

Retentie

Of de student doorstudeert na het eerste studiejaar

Ja, Nee

SES_Arbeid

De sociaaleconomische status score op arbeid van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding

-1 - 1

SES_Totaal

De sociaaleconomische status score totaal van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding

-1 - 1

SES_Welvaart

De sociaaleconomische status score op welvaart van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding

-1 - 1

Studiekeuzeprofiel

Het studiekeuzeprofiel in het voortgezet onderwijs (havo of vwo) of profiel in het mbo

AHO, ALG, CERT, CM, EA, EM, EM&CM, HB, HO, ICT, MedV, NG, NT, NT&NG, Onbekend, TP, TR, TSL, VNL, VS, ZW

Vooropleiding

De toelaatgevende vooropleiding

BD, CD, HAVO, HO, MBO, Overig, VWO

1.1.3 Toelichting op de analyse

We toetsen in deze analyse Retentie na 1 jaar, voortaan Retentie genoemd.

Retentie is gedefinieerd als ingeschreven staan in dezelfde opleiding in een aansluitend collegejaar. Een wisseling van opleidingsvorm binnen de opleiding, bijvoorbeeld van voltijd in jaar 1 naar duaal in jaar 2, geldt ook als retentie. Uitval is het tegenovergestelde van retentie: niet ingeschreven staan in dezelfde opleiding in een aansluitend collegejaar. Een wisseling van opleidingsvorm binnen de opleiding, bijvoorbeeld van voltijd in jaar 1 naar duaal in jaar 2, geldt niet als uitval.

1.2 Voorbereidingen

1.2.1 Laad de data

We laden een subset in van historische data specifiek voor:

Opleiding: ITD | B Communication and Multimedia Design (Synth) (CMD), voltijd, eerstejaars - Retentie na 1 jaar

Toon code
## Read the data for this study progamme
if(params$use_synthetic_data) {
  dfOpleiding_inschrijvingen_base <- Get_Studyprogram_Enrollments_Synthetic(
    studytrack = params$opleiding,
    studyform = toupper(params$opleidingsvorm_afkorting)
  ) |> 
    mutate(
      INS_Student_UUID_opleiding_vorm = paste(ID, INS_Opleiding, INS_Opleidingsvorm, sep = "_"),
      INS_Opleidingsnaam_huidig = paste(INS_Opleidingsnaam_huidig, "(Synth.)", sep = " ")
    )
} else {
  dfOpleiding_inschrijvingen_base <- get_lta_studyprogram_enrollments_pin(
    board = "HHs/Inschrijvingen",
    faculty = params$faculteit,
    studyprogram = params$opleidingsnaam,
    studytrack = params$opleiding,
    studyform = toupper(params$opleidingsvorm),
    range = "eerstejaars")
}

## Rearrange the levels
Set_Levels(dfOpleiding_inschrijvingen_base)

dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |>  
  
  ## Create a simple success variable
  Mutate_Retentie(sSucces_model) |>
  
  ## Convert the success variable into a factor
  mutate(SUC_Retentie = as.factor(SUC_Retentie)) |> 

  ## Special possibly based on the propaedeutic diploma
  # Filter_Propedeusediploma(sPropedeusediploma) |>

  ## Make the Dual Study variable a Yes/No variable
  mutate(INS_Dubbele_studie = ifelse(INS_Aantal_inschrijvingen > 1, "Ja", "Nee")) |>  

  ## Remove INS_Aantal_inschrijvingen
  select(-INS_Aantal_inschrijvingen) |> 

  ## Adjust levels for some variables
  Mutate_Levels(
  c(
    "VOP_Studiekeuzeprofiel_LTA_afkorting",
    "INS_Aansluiting_LTA",
    "VOP_Toelaatgevende_vooropleiding_soort"
  ),
    list(lLevels_skp, lLevels_vop, lLevels_vop)
  )
  
## B Huidtherapie: Filter on only students with a grade number (selection)
if(opleiding == "HDT") {
  dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |> 
    filter(!is.na(RNK_Rangnummer)) 
} 

1.2.2 Selecteer en inspecteer de data

We selecteren eerst de relevante variabelen. We verwijderen daarbij variabelen die maar 1 waarde hebben. We inspecteren de variabelen in een samenvatting in relatie tot retentie en corrigeren daarbij voor multiple testing; de gecorrigeerde significantie-waarden staan vermeld als q-value. Daarnaast inspecteren we de kwaliteit van de data op missende waarden.

Toon code
lSelect <- c(
    "INS_Student_UUID_opleiding_vorm",
    "CBS_APCG_tf",
    "DEM_Geslacht",
    "DEM_Leeftijd_1_oktober",
    "GIS_Tijd_fiets_OV",
    "INS_Collegejaar",
    "INS_Dagen_tussen_aanmelding_en_1_september",
    "INS_Dubbele_studie",
    "INS_Aansluiting_LTA",
    "SES_Deelscore_arbeid",
    "SES_Deelscore_welvaart",
    "SES_Totaalscore",
    "SUC_Retentie",
    "VOP_Gemiddeld_cijfer_cijferlijst",
    "VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO",
    "VOP_Cijfer_CE1_nederlands",
    "VOP_Cijfer_CE1_engels",
    "VOP_Cijfer_CE_proxy_wiskunde",
    "VOP_Cijfer_CE1_natuurkunde",
    "VOP_Studiekeuzeprofiel_LTA_afkorting",
    "VOP_Toelaatgevende_vooropleiding_soort"
  )

## B Huidtherapie: add the variable RNK_Rangnummer
if(opleiding == "HDT") {
  lSelect <- c(lSelect, "RNK_Rangnummer")
}

## Create a subset
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen_base |>
  
  ## Select the relevant variables
  select_at(lSelect) |>
  
  ## Rename variables for more readable names
  rename(
    ID                    = INS_Student_UUID_opleiding_vorm,
    Geslacht              = DEM_Geslacht,
    Leeftijd              = DEM_Leeftijd_1_oktober,
    Reistijd              = GIS_Tijd_fiets_OV,
    Dubbele_studie        = INS_Dubbele_studie,
    Collegejaar           = INS_Collegejaar,
    Aanmelding            = INS_Dagen_tussen_aanmelding_en_1_september,
    Aansluiting           = INS_Aansluiting_LTA,
    APCG                  = CBS_APCG_tf,
    SES_Arbeid            = SES_Deelscore_arbeid,
    SES_Welvaart          = SES_Deelscore_welvaart,
    SES_Totaal            = SES_Totaalscore,          
    Retentie              = SUC_Retentie,
    Cijfer_SE_VO          = VOP_Gemiddeld_cijfer_cijferlijst,
    Cijfer_CE_VO          = VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO,
    Cijfer_CE_Nederlands  = VOP_Cijfer_CE1_nederlands,
    Cijfer_CE_Engels      = VOP_Cijfer_CE1_engels,
    Cijfer_CE_Wiskunde    = VOP_Cijfer_CE_proxy_wiskunde,
    Cijfer_CE_Natuurkunde = VOP_Cijfer_CE1_natuurkunde,
    Studiekeuzeprofiel    = VOP_Studiekeuzeprofiel_LTA_afkorting,
    Vooropleiding         = VOP_Toelaatgevende_vooropleiding_soort
  ) |> 
  
  ## Adjust CBS_APCG_tf to a factor
  mutate(APCG = case_when(APCG == TRUE ~ "Ja",
                          APCG == FALSE ~ "Nee",
                          .default = "Onbekend")) |>

  ## Indicate where missing numbers are in VO
  Mutate_Cijfers_VO() |>
  
  ## Remove variables, where there is only 1 value
  select(where(~ n_distinct(.) > 1)) |>
  
  arrange(Collegejaar, ID)

## B Huidtherapie: rename the variable RNK_Rangnummer
if(opleiding == "HDT") {
  dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
    rename(Rangnummer = RNK_Rangnummer)
} 

dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
 ltabase::sort_distinct()

## Remove the base dataset
## rm(dfOpleiding_inschrijvingen_base)
Tabel 1.2: Variabelen in relatie tot de uitkomstmaat: params$succes

Retentie

Variabele

Ja, N=1002 (62%)1

Nee, N=611 (38%)1

p-value2

q-value3

Totaal, N = 16131

Aanmelding

135,80 (62,44)

112,36 (68,57)

<0,001

<0,001***

126,92 (65,80)

Aansluiting

<0,001

0,012***

2e Studie

6 (40%)

9 (60%)

15 (100%)

Direct

455 (60%)

301 (40%)

756 (100%)

Na CD

14 (74%)

5 (26%)

19 (100%)

Onbekend

0 (NA%)

0 (NA%)

0 (NA%)

Overig

0 (NA%)

0 (NA%)

0 (NA%)

Switch extern

254 (57%)

191 (43%)

445 (100%)

Switch intern

148 (73%)

54 (27%)

202 (100%)

Tussenjaar

125 (71%)

51 (29%)

176 (100%)

APCG

<0,001

0,013***

Ja

224 (63%)

133 (37%)

357 (100%)

Nee

721 (64%)

411 (36%)

1.132 (100%)

Onbekend

57 (46%)

67 (54%)

124 (100%)

Cijfer_CE_Engels

6,93 (1,18)

7,07 (1,10)

0,16

>0,99

6,98 (1,15)

Cijfer_CE_Engels_missing

0,70

>0,99

Ja

533 (62%)

331 (38%)

864 (100%)

Nee

469 (63%)

280 (37%)

749 (100%)

Cijfer_CE_Natuurkunde

6,18 (0,96)

6,31 (0,95)

0,21

>0,99

6,23 (0,95)

Cijfer_CE_Natuurkunde_missing

0,68

>0,99

Ja

867 (62%)

533 (38%)

1.400 (100%)

Nee

135 (63%)

78 (37%)

213 (100%)

Cijfer_CE_Nederlands

5,94 (0,85)

5,96 (0,91)

0,86

>0,99

5,95 (0,87)

Cijfer_CE_Nederlands_missing

0,61

>0,99

Ja

533 (62%)

333 (38%)

866 (100%)

Nee

469 (63%)

278 (37%)

747 (100%)

Cijfer_CE_VO

6,46 (0,37)

6,35 (0,34)

<0,001

<0,001***

6,42 (0,36)

Cijfer_CE_VO_missing

0,59

>0,99

Ja

488 (61%)

306 (39%)

794 (100%)

Nee

514 (63%)

305 (37%)

819 (100%)

Cijfer_CE_Wiskunde

6,35 (1,12)

6,28 (1,07)

0,65

>0,99

6,33 (1,10)

Cijfer_CE_Wiskunde_missing

0,64

>0,99

Ja

544 (62%)

339 (38%)

883 (100%)

Nee

458 (63%)

272 (37%)

730 (100%)

Cijfer_SE_VO

6,46 (0,38)

6,35 (0,36)

<0,001

<0,001***

6,42 (0,38)

Cijfer_SE_VO_missing

0,59

>0,99

Ja

501 (61%)

314 (39%)

815 (100%)

Nee

501 (63%)

297 (37%)

798 (100%)

Dubbele_studie

<0,001

<0,001***

Ja

12 (26%)

34 (74%)

46 (100%)

Nee

990 (63%)

577 (37%)

1.567 (100%)

Geslacht

<0,001

<0,001***

M

521 (57%)

394 (43%)

915 (100%)

V

481 (69%)

217 (31%)

698 (100%)

Leeftijd

20,01 (2,26)

20,19 (2,61)

0,48

>0,99

20,08 (2,40)

Reistijd

37,01 (18,38)

38,84 (22,79)

0,79

>0,99

37,70 (20,17)

SES_Arbeid

0,01 (0,08)

0,01 (0,08)

0,57

>0,99

0,01 (0,08)

SES_Totaal

0,02 (0,28)

0,01 (0,27)

0,45

>0,99

0,02 (0,28)

SES_Welvaart

0,01 (0,13)

0,01 (0,13)

0,36

>0,99

0,01 (0,13)

Studiekeuzeprofiel

0,054

>0,99

AHO

2 (67%)

1 (33%)

3 (100%)

ALG

0 (0%)

1 (100%)

1 (100%)

CERT

1 (100%)

0 (0%)

1 (100%)

CM

78 (72%)

30 (28%)

108 (100%)

EA

36 (67%)

18 (33%)

54 (100%)

EM

118 (56%)

93 (44%)

211 (100%)

EM&CM

152 (70%)

66 (30%)

218 (100%)

HB

8 (62%)

5 (38%)

13 (100%)

HO

27 (57%)

20 (43%)

47 (100%)

ICT

82 (64%)

47 (36%)

129 (100%)

MedV

132 (65%)

72 (35%)

204 (100%)

NG

115 (63%)

68 (37%)

183 (100%)

NT

56 (60%)

38 (40%)

94 (100%)

NT&NG

62 (60%)

41 (40%)

103 (100%)

Onbekend

90 (52%)

83 (48%)

173 (100%)

TP

2 (40%)

3 (60%)

5 (100%)

TR

5 (56%)

4 (44%)

9 (100%)

TSL

4 (80%)

1 (20%)

5 (100%)

VNL

8 (80%)

2 (20%)

10 (100%)

VS

1 (33%)

2 (67%)

3 (100%)

ZW

23 (59%)

16 (41%)

39 (100%)

Vooropleiding

0,010

0,25*

BD

42 (46%)

50 (54%)

92 (100%)

CD

21 (70%)

9 (30%)

30 (100%)

HAVO

548 (64%)

312 (36%)

860 (100%)

HO

27 (53%)

24 (47%)

51 (100%)

MBO

330 (63%)

192 (37%)

522 (100%)

Overig

0 (NA%)

0 (NA%)

0 (NA%)

VWO

34 (59%)

24 (41%)

58 (100%)

1Mean (SD); n (%)

2Wilcoxon rank sum test; Fisher's Exact Test for Count Data with simulated p-value
(based on 2000 replicates); Pearson's Chi-squared test

3*p<0.05; **p<0.01; ***p<0.001

Toon code
## Load dlookr (temporary)
suppressMessages(library(dlookr))

## Show a summary of the data, sorted by missing values
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(missing_percent, 2)) |>
  arrange(desc(missing_percent)) |>
  knitr::kable(col.names = c("Variabelen",
                           "Type",
                           "# Missende waarden",
                           "% Missende waarden",
                           "# Unieke waarden",
                           "Ratio unieke waarden"))
## Detach dlookr
detach("package:dlookr", unload = TRUE)
Tabel 1.3: Kwaliteit van de data voor bewerkingen (gesorteerd op missende waarden)
Variabelen Type # Missende waarden % Missende waarden # Unieke waarden Ratio unieke waarden
Cijfer_CE_Natuurkunde numeric 1400 86.79 42 86.79
Cijfer_CE_Wiskunde numeric 883 54.74 56 54.74
Cijfer_CE_Nederlands numeric 866 53.69 45 53.69
Cijfer_CE_Engels numeric 864 53.56 57 53.56
Cijfer_SE_VO numeric 815 50.53 23 50.53
Cijfer_CE_VO numeric 794 49.23 23 49.23
Studiekeuzeprofiel factor 173 10.73 21 10.73
SES_Welvaart numeric 126 7.81 355 7.81
SES_Arbeid numeric 125 7.75 261 7.75
SES_Totaal numeric 125 7.75 459 7.75
Reistijd numeric 20 1.24 362 1.24
Aanmelding numeric 0 0.00 266 0.00
Aansluiting factor 0 0.00 6 0.00
APCG character 0 0.00 3 0.00
Cijfer_CE_Engels_missing character 0 0.00 2 0.00
Cijfer_CE_Natuurkunde_missing character 0 0.00 2 0.00
Cijfer_CE_Nederlands_missing character 0 0.00 2 0.00
Cijfer_CE_VO_missing character 0 0.00 2 0.00
Cijfer_CE_Wiskunde_missing character 0 0.00 2 0.00
Cijfer_SE_VO_missing character 0 0.00 2 0.00
Collegejaar numeric 0 0.00 11 0.00
Dubbele_studie character 0 0.00 2 0.00
Geslacht factor 0 0.00 2 0.00
ID character 0 0.00 1613 0.00
Leeftijd numeric 0 0.00 17 0.00
Retentie factor 0 0.00 2 0.00
Vooropleiding factor 0 0.00 6 0.00

1.2.3 Bewerk de data

  • Uit de eerste diagnose blijkt dat niet alle variabelen goed genoeg zijn voor het bouwen van een prognosemodel: er zijn missende waarden en niet alle veldtypes zijn geschikt.
  • Om bias te voorkomen verwijderen we geen rijen met missende waarden, maar vullen die op (imputatie). We bewerken de data zo dat alle missende waarden worden opgevuld: bij numerieke waarden met het gemiddelde en bij categorische variabelen met ‘Onbekend’.
  • We passen het type van sommige variabelen aan, zodat ze in het model gebruikt kunnen worden: tekstvelden zetten we om naar factor (een categorische variabele); logische variabelen (Ja/Nee) zetten we om naar een numerieke variabele (1/0).
  • De uitkomstvariabele, Retentie, leiden we af van de variabele SUC_Uitval_aantal_jaar_LTA. Als de waarde daar 1 is, is de student na 1 jaar uitgevallen, 2 na 2 jaar, etc. Zolang de waarde daar 0 is, is de student niet uitgevallen.
  • Een fictief studentnummer (INS_Student_UUID_opleiding_vorm) gebruiken we, zodat we - als er afwijkende resultaten zijn - de dataset gericht kunnen onderzoeken als dat nodig is.
Toon code
## Edit the data
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
  
  ## Imputate all numeric variables with the mean
  mutate(across(where(is.numeric), ~ ifelse(
    is.na(.x),
    mean(.x, na.rm = T),
    .x
  )) ) |>
  
  ## Convert character variables to factor
  mutate(across(where(is.character), as.factor)) |> 
  
  ## Convert logical variables to 0 or 1
  mutate(across(where(is.logical), as.integer)) |>
  
  ## Fill in factors missing values with “Unknown”
  mutate(across(where(is.factor), ~ suppressWarnings(
    fct_explicit_na(.x, na_level = "Onbekend")
  ))) |> 
  
  ## Rearrange the columns so that Retentie is in front
  select(Retentie, everything()) 

## View the data
## glimpse(dfOpleiding_inschrijvingen) 

## Load dlookr (temporary)
suppressMessages(library(dlookr))

## Diagnose the data
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(unique_rate, 2)) |>
  knitr::kable(col.names = c("Variabelen",
                           "Type",
                           "# Missende waarden",
                           "% Missende waarden",
                           "# Unieke waarden",
                           "Ratio unieke waarden"))
## Detach dlookr
detach("package:dlookr", unload = TRUE)
Tabel 1.4: Kwaliteit van de data na bewerkingen (gesorteerd op missende waarden)
Variabelen Type # Missende waarden % Missende waarden # Unieke waarden Ratio unieke waarden
Retentie factor 0 0 2 0.00
Aanmelding numeric 0 0 266 0.16
Aansluiting factor 0 0 6 0.00
APCG factor 0 0 3 0.00
Cijfer_CE_Engels numeric 0 0 57 0.04
Cijfer_CE_Engels_missing factor 0 0 2 0.00
Cijfer_CE_Natuurkunde numeric 0 0 42 0.03
Cijfer_CE_Natuurkunde_missing factor 0 0 2 0.00
Cijfer_CE_Nederlands numeric 0 0 45 0.03
Cijfer_CE_Nederlands_missing factor 0 0 2 0.00
Cijfer_CE_VO numeric 0 0 23 0.01
Cijfer_CE_VO_missing factor 0 0 2 0.00
Cijfer_CE_Wiskunde numeric 0 0 56 0.03
Cijfer_CE_Wiskunde_missing factor 0 0 2 0.00
Cijfer_SE_VO numeric 0 0 23 0.01
Cijfer_SE_VO_missing factor 0 0 2 0.00
Collegejaar numeric 0 0 11 0.01
Dubbele_studie factor 0 0 2 0.00
Geslacht factor 0 0 2 0.00
ID factor 0 0 1613 1.00
Leeftijd numeric 0 0 17 0.01
Reistijd numeric 0 0 362 0.22
SES_Arbeid numeric 0 0 261 0.16
SES_Totaal numeric 0 0 459 0.28
SES_Welvaart numeric 0 0 355 0.22
Studiekeuzeprofiel factor 0 0 21 0.01
Vooropleiding factor 0 0 6 0.00

1.2.4 Bekijk de onderlinge correlaties

Het is verstandig om voorafgaand aan het bouwen van een model te kijken naar de onderlinge correlaties tussen numerieke variabelen. Dit geeft inzicht in de data en kan helpen bij het maken van keuzes voor het model of de duiding van de uitkomsten.

Toon code
## Create a plot of the intercorrelations in numerical variables
## Remove columns with a standard deviation of 0
dfOpleiding_inschrijvingen |> 
  select(-Collegejaar) |>
  select(where(is.numeric)) |> 
  select_if(~ sd(.) != 0) |>
  cor() |> 
  corrplot::corrplot(
    order = 'hclust', 
    addrect = 4,
    method = "number",  
    tl.cex = 0.8,       
    tl.col = "black",
    diag = FALSE)
Figuur 1.1: Correlatiematrix

1.2.5 Bouw de trainingset, validatieset en testset

  • De data is nu geschikt om een prognosemodel mee te bouwen.
  • Om het model te bouwen, testen en valideren, splitsen we de data in drie delen van 60%, 20% en 20%. We doen dit op zo’n manier, dat elk deel ongeveer een gelijk aantal studenten bevat dat doorstudeert (dus niet uitvalt).
  • We trainen het model op basis van 60% en valideren de modellen tijdens het trainen op de overige 20% (de validatieset).
  • De verdeling van de training- en validatieset muteren we 10x (10 folds) om te voorkomen dat het model te veel leert van de trainingset en daardoor slecht presteert op de validatieset.
  • Als het model klaar is, testen we het op de 20% studenten uit de testset. De testset blijft dus de gehele tijd ongemoeid, zodat we overfitting - een te goed model op bekende data, maar slechte presetaties (performance) op onbekende data - voorkomen.
  • Een willekeurig, maar vaststaand seed-getal voorkomt dat we bij elke run van het model c.q. deze code een net iets andere uitkomst krijgen.
Toon code
knitr::include_graphics(here::here("R/images", "voorspelmodel-dataset-lta-hhs.png"))
Figuur 1.2: Splitsing van de dataset in trainingset, validatieset en testset
Toon code
set.seed(0821)

## Split the data into 3 parts: 60%, 20% and 20%
splits      <- initial_validation_split(dfOpleiding_inschrijvingen,
                                        strata = Retentie,
                                        prop = c(0.6, 0.2))

## Create three sets: a training set, a test set and a validation set
dfRetentie_train      <- training(splits)
dfRetentie_test       <- testing(splits)
dfRetentie_validation <- validation_set(splits)

## Create a resample set based on 10 folds (default)
dfRetentie_resamples  <- vfold_cv(dfRetentie_train, strata = Retentie)
Tabel 1.5: Verhouding van de uitkomstvariabele in de training- en testset
Naam Retentie Aantal Proportie
Trainingset FALSE 366 37.8%
Trainingset TRUE 601 62.2%
Testset FALSE 123 38.0%
Testset TRUE 201 62.0%

1.3 Model I: Logistische Regressie

  • Het eerste model is een logistische regressie met penalized likelihood; we gebruiken de glmnet engine voor het bouwen van het model. Penalized likelihood is een techniek die helpt bij het voorkomen van overfitting: het voegt voor elke extra variabele een strafterm toe om eenvoudige modellen te belonen. Glmnet is een veelgebruikt package voor het bouwen van logistische regressiemodellen.
  • We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. De ROC-curve (Receiver Operating Characteristic) is een grafiek die de prestaties van een classificatiemodel afbeeldt door de verhouding tussen de true positives (sensitiviteit) en de false positives (aspecficiteit = 1-specificiteit) te plotten bij verschillende drempelwaarden. De oppervlakte onder deze curve, bekend als de AUC (Area Under the Curve), kwantificeert het onderscheidingsvermogen van het model; een AUC van 1 duidt op een perfect onderscheidend vermogen, terwijl een AUC van 0,5 wijst op een model zonder onderscheidend vermogen.

1.3.1 Maak het model

Eerst bouwen we het model.

## Build the model: logistic regression
lr_mod <- 
  logistic_reg(penalty = tune(), mixture = 1) |> 
  set_engine("glmnet")

1.3.2 Maak de recipe

Vervolgens zetten we meerdere stappen in een ‘recipe’:

  • We definiëren de student-ID als ID variabele. Daarmee krijgt deze variabele de rol van uniek rij-kenmerk.
  • We verwijderen vervolgens de oorspronkelijke student-ID en het collegejaar uit de data, omdat deze verder niet gebruikt moeten worden in het model.
  • We converteren factoren naar dummy variabelen: voor elke categorie wordt er een nieuwe logische variabele (Ja/Nee) aangemaakt.
  • We verwijderen variabelen die geen waarde toevoegen: variabelen met uitsluitend nullen.
  • We normaliseren numerieke variabelen om ze met elkaar te kunnen vergelijken door ze te centreren en schalen.
  • Sterk gecorreleerde waarden verwijderen we nu niet, omdat we later in de analyse de eventuele samenhang met andere variabelen in een prognosemodel nog willen kunnen visualiseren.
## Build the recipe: logistic regression
lr_recipe <- 
  recipe(Retentie ~ ., data = dfRetentie_train) |>  
  update_role(ID, new_role = "ID") |>           ## Set the student ID as an ID variable
  step_rm(ID, Collegejaar) |>                   ## Remove ID and college year from the model
  step_dummy(all_nominal_predictors()) |>       ## Create dummy variables from categorical variables
  step_zv(all_predictors()) |>                  ## Remove zero values
  step_normalize(all_numeric_predictors())      ## Center and scale numeric variables
Toon code
## Show the recipe
tidy(lr_recipe) |> 
  knitr::kable(col.names = c("Nummer", 
                             "Operatie", 
                             "Type",
                             "Getraind",
                             "Sla over",
                             "ID"))
Tabel 1.6: Recipesteps voor logistische regressie
Nummer Operatie Type Getraind Sla over ID
1 step rm FALSE FALSE rm_KV8Lp
2 step dummy FALSE FALSE dummy_bUMsj
3 step zv FALSE FALSE zv_beaLb
4 step normalize FALSE FALSE normalize_3WTmw

De variabelen die nu nog overblijven zijn:

Tabel 1.7: Resterende variabelen voor logistische regressie na bewerkingen
Aanmelding APCG_Nee Studiekeuzeprofiel_MedV
Cijfer_CE_Engels APCG_Onbekend Studiekeuzeprofiel_NG
Cijfer_CE_Natuurkunde Cijfer_CE_Engels_missing_Nee Studiekeuzeprofiel_NT
Cijfer_CE_Nederlands Cijfer_CE_Natuurkunde_missing_Nee Studiekeuzeprofiel_NT.NG
Cijfer_CE_VO Cijfer_CE_Nederlands_missing_Nee Studiekeuzeprofiel_Onbekend
Cijfer_CE_Wiskunde Cijfer_CE_VO_missing_Nee Studiekeuzeprofiel_TP
Cijfer_SE_VO Cijfer_CE_Wiskunde_missing_Nee Studiekeuzeprofiel_TR
Leeftijd Cijfer_SE_VO_missing_Nee Studiekeuzeprofiel_TSL
Reistijd Dubbele_studie_Nee Studiekeuzeprofiel_VNL
SES_Arbeid Geslacht_V Studiekeuzeprofiel_VS
SES_Totaal Studiekeuzeprofiel_ALG Studiekeuzeprofiel_ZW
SES_Welvaart Studiekeuzeprofiel_CM Vooropleiding_CD
Retentie Studiekeuzeprofiel_EA Vooropleiding_HAVO
Aansluiting_Direct Studiekeuzeprofiel_EM Vooropleiding_HO
Aansluiting_Na.CD Studiekeuzeprofiel_EM.CM Vooropleiding_MBO
Aansluiting_Switch.extern Studiekeuzeprofiel_HB Vooropleiding_VWO
Aansluiting_Switch.intern Studiekeuzeprofiel_HO
Aansluiting_Tussenjaar Studiekeuzeprofiel_ICT

1.3.3 Maak de workflow

Voor de uitvoering bouwen we een workflow. Daaraan voegen we het model en de bewerkingen in de recipe toe.

## Create the workflow: logistic regression
lr_workflow <- 
  workflow() |>         ## Create a workflow
  add_model(lr_mod) |>  ## Add the model
  add_recipe(lr_recipe) ## Add the recipe

## Show workflow
lr_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
4 Recipe Steps

• step_rm()
• step_dummy()
• step_zv()
• step_normalize()

── Model ───────────────────────────────────────────────────────────────────────
Logistic Regression Model Specification (classification)

Main Arguments:
  penalty = tune()
  mixture = 1

Computational engine: glmnet 

1.3.4 Tune en train het model

Het model moet getuned worden. Dit houdt in dat we de beste parameters voor het model moeten vinden. We maken een grid met verschillende penalty waarden. Daarmee kunnen we vervolgens het beste model selecteren met de hoogste ROC/AUC. We plotten de resultaten van de tuning, zodat we hieruit het beste model kunnen kiezen.

## Create a grid: logistic regression
lr_reg_grid <- tibble(penalty = 10 ^ seq(-4, -1, length.out = 30))

## Train and tune the model: logistic regression
lr_res <- 
  lr_workflow |> 
  tune_grid(dfRetentie_validation,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
Toon code
## Plot the results + a red vertical line for the max AUC
lr_plot <- 
  lr_res |> 
  collect_metrics() |> 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  
  ## Make the scale of the x-axis logarithmic
  scale_x_log10(labels = scales::label_number()) +
    theme(
      axis.title.x = element_text(margin = margin(t = 20))
    ) +
  
  # Define the title, subtitle and caption
  labs(
    caption = sCaption,
    x = "Area under the ROC Curve",
    y = "Penalty"
  )
  
  ## Add LTA elements
  lr_plot <- Add_LTA_Theme_Elements(lr_plot, title_subtitle = FALSE)
  
# Find the penalty value with the max AUC
max_auc_penalty <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)

# Add the red vertical line to lr_plot
lr_plot_plus <- lr_plot + 
  geom_vline(xintercept = max_auc_penalty, color = "red")

# Find a mean for the max AUC that is higher
max_auc_mean <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)

## Print the final plot
lr_plot_plus
Figuur 1.3: Tuning resultaten logistische regressie

1.3.5 Kies het beste model

De prestaties van een model gevisualiseerd met behulp van een ROC curve. De sensitiviteit (True Positive Rate) en specificiteit (True Negative Rate) worden hierin tegenover elkaar uitgezet. De Area under the ROC Curve (AUC/ROC) geeft de prestaties van het model weer. Het model scoort beter naarmate de AUC/ROC dichter bij de 1 ligt, de linker bovenhoek. De linker bovenhoek houdt in dat alle prognoses exact overeenstemmen met de werkelijkheid. Een AUC/ROC van 0,5 betekent dat het model niet beter presteert dan een willekeurige voorspelling.

We gebruiken modellen met een zo hoog mogelijke Area under the ROC Curve (AUC/ROC) en een zo laag mogelijke penalty. Zo kunnen we uit de resultaten het beste model kiezen en visualiseren.

## Show the best model
top_models <-
  lr_res |> 
  show_best(metric = "roc_auc", n = 10) |> 
  mutate(mean = round(mean, 6)) |>
  arrange(penalty) 
Toon code
top_models|> 
  knitr::kable(col.names = c("Penalty", 
                             "Metriek", 
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))
Tabel 1.8: Model performance voor logistische regressie
Penalty Metriek Estimator Gemiddelde Aantal SE Configuratie
0.0017433 roc_auc binary 0.668279 1 NA Preprocessor1_Model13
0.0022122 roc_auc binary 0.670205 1 NA Preprocessor1_Model14
0.0028072 roc_auc binary 0.670820 1 NA Preprocessor1_Model15
0.0035622 roc_auc binary 0.672172 1 NA Preprocessor1_Model16
0.0045204 roc_auc binary 0.673934 1 NA Preprocessor1_Model17
0.0057362 roc_auc binary 0.675369 1 NA Preprocessor1_Model18
0.0072790 roc_auc binary 0.676967 1 NA Preprocessor1_Model19
0.0092367 roc_auc binary 0.678934 1 NA Preprocessor1_Model20
0.0117210 roc_auc binary 0.677336 1 NA Preprocessor1_Model21
0.0148735 roc_auc binary 0.669262 1 NA Preprocessor1_Model22
## Select the best model: logistic regression
lr_best <- 
  lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |>
  slice(1) 
Toon code
lr_best|> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable(col.names = c("Penalty", 
                             "Metriek", 
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))
Tabel 1.9: Hoogste model performance voor logistische regressie
Penalty Metriek Estimator Gemiddelde Aantal SE Configuratie
0.0092367 roc_auc binary 0.678934 1 NA Preprocessor1_Model20
## Collect the predictions and evaluate the model (AUC/ROC): logistic regression
lr_auc <- 
  lr_res |> 
  collect_predictions(parameters = lr_best) |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Logistisch Regressie")
Toon code
## Plot the ROC curve
Get_ROC_Plot(lr_auc, position = 1)
Figuur 1.4: ROC curve voor logistische regressie
Toon code
## Determine the AUC of the best model
lr_auc_highest   <-
  lr_res |>
  collect_predictions(parameters = lr_best) |> 
  roc_auc(Retentie, .pred_FALSE)

## Add model name and AUC dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Logistic Regression", auc = lr_auc_highest$.estimate)

1.4 Model II: Tree-based ensemble

  • Het tweede model is een random forest: een ensemble van beslisbomen (decision trees). Het is een krachtig model dat goed om kan gaan met complexe data en veel variabelen.
  • We gebruiken de ranger engine voor het bouwen van het model.

1.4.1 Bepaal het aantal PC-cores

Omdat een random forest model veel berekeningen vereist, willen we daarvoor alle computerkracht gebruiken die beschikbaar is. Het aantal CPU’s (cores), wat verschilt per computer, bepaalt hoe snel het model getraind kan worden. We bepalen het aantal cores en gebruiken dat bij het bouwen van het model.

Toon code
## Determine the number of cores
cores <- parallel::detectCores()

1.4.2 Maak het model

We bouwen eerst het model. We gebruiken de rand_forest functie om het model te bouwen. We tunen de mtry en min_n parameters. De mtry parameter bepaalt het aantal variabelen dat per boom wordt gebruikt. De min_n parameter bepaalt het minimum aantal observaties dat in een blad van de boom moet zitten. De functie tune() is hier nog een placeholder om de beste waarden voor deze parameters - die we later bepalen - in te kunnen stellen. We gebruiken 1.000 bomen c.q. versies van het model.

## Build the model: random forest

rf_mod <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) |> 
  set_engine("ranger", num.threads = cores) |> 
  set_mode("classification")

1.4.3 Maak de recipe

We maken een recipe voor het random forest model. We verwijderen de student ID en het collegejaar uit de data, omdat deze niet moet worden gebruikt in het model. Overige stappen zijn bij een random forest minder relevant in tegenstelling tot een regressiemodel.

## Create the recipe: random forest
rf_recipe <- 
  recipe(Retentie ~ ., data = dfRetentie_train) |> 
  step_rm(ID, Collegejaar)                      ## Verwijder ID en Collegejaar uit het model
Toon code
## Show the recipe
tidy(rf_recipe) |> 
  knitr::kable(col.names = c("Nummer", 
                             "Operatie", 
                             "Type",
                             "Getraind",
                             "Sla over",
                             "ID"))
Tabel 1.10: Recipesteps voor random forest
Nummer Operatie Type Getraind Sla over ID
1 step rm FALSE FALSE rm_HCcZp

1.4.4 Maak de workflow

We voegen het model en de recipe toe aan de workflow voor dit model.

## Create the workflow: random forest
rf_workflow <- 
  workflow() |> 
  add_model(rf_mod) |> 
  add_recipe(rf_recipe)

## Show workflow
rf_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ────────────────────────────────────────────────────────────────
1 Recipe Step

• step_rm()

── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger 

1.4.5 Tune en train het model

We trainen en tunen het model in de workflow. We maken een grid met verschillende waarden voor de parameters mtry en min_n. We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. Met de resultaten van de tuning kiezen we het beste model.

## Show the parameters that can be tuned
rf_mod
Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger 
## Extract the parameters being tuned
extract_parameter_set_dials(rf_mod)
Collection of 2 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[?]
      min_n min_n nparam[+]

Model parameters needing finalization:
   # Randomly Selected Predictors ('mtry')

See `?dials::finalize` or `?dials::update.parameters` for more information.
## Determine the seed
set.seed(2904)

## Build the grid: random forest
rf_res <- 
  rf_workflow |> 
  tune_grid(dfRetentie_validation,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
i Creating pre-processing data to finalize unknown parameter: mtry

1.4.6 Kies het beste model

We evalueren de beste modellen en maken een ROC curve om de performance van het model te visualiseren. Vervolgens vergelijken we de prestaties van de modellen en kiezen daaruit het beste model.

## Show the best models
rf_res |> 
  show_best(metric = "roc_auc", n = 15) |> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable(col.names = c("Mtry", 
                             "Min. aantal", 
                             "Metriek",
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))
Tabel 1.11: Model performance voor random forest
Mtry Min. aantal Metriek Estimator Gemiddelde Aantal SE Configuratie
2 6 roc_auc binary 0.700041 1 NA Preprocessor1_Model22
11 29 roc_auc binary 0.699836 1 NA Preprocessor1_Model02
4 26 roc_auc binary 0.699385 1 NA Preprocessor1_Model23
5 36 roc_auc binary 0.699344 1 NA Preprocessor1_Model09
7 35 roc_auc binary 0.698074 1 NA Preprocessor1_Model14
8 31 roc_auc binary 0.698033 1 NA Preprocessor1_Model24
11 21 roc_auc binary 0.697951 1 NA Preprocessor1_Model15
2 32 roc_auc binary 0.696434 1 NA Preprocessor1_Model03
3 8 roc_auc binary 0.695492 1 NA Preprocessor1_Model08
6 23 roc_auc binary 0.695246 1 NA Preprocessor1_Model17
13 25 roc_auc binary 0.694262 1 NA Preprocessor1_Model25
22 38 roc_auc binary 0.692295 1 NA Preprocessor1_Model01
16 34 roc_auc binary 0.691926 1 NA Preprocessor1_Model07
18 40 roc_auc binary 0.690369 1 NA Preprocessor1_Model16
15 20 roc_auc binary 0.690287 1 NA Preprocessor1_Model20
Toon code
## Plot the results
autoplot <- autoplot(rf_res) +
  theme_minimal() +
  labs(
    y = "roc/auc",
    caption = sCaption
  )
  
  ## Add LTA elements
  autoplot <- Add_LTA_Theme_Elements(autoplot, title_subtitle = FALSE)
  
  print(autoplot)
Figuur 1.5: Model performance random forest
## Select the best model
rf_best <- 
  rf_res |> 
  select_best(metric = "roc_auc")
Toon code
rf_best|> 
  knitr::kable(col.names = c("Mtry", 
                             "Min. aantal", 
                             "Configuratie"))
Tabel 1.12: Hoogste model performance voor random forest
Mtry Min. aantal Configuratie
2 6 Preprocessor1_Model22
Toon code
## Collect the predictions
rf_res |> 
  collect_predictions() |> 
  head(10) |>
  mutate(.pred_FALSE = scales::percent(.pred_FALSE, accuracy = 0.1),
         .pred_TRUE = scales::percent(.pred_TRUE, accuracy = 0.1)) |>
  knitr::kable(col.names = c("% Voorsp. FALSE", 
                             "% Voorsp. TRUE", 
                             "ID",
                             "Rij",
                             "Mtry", 
                             "Min. aantal", 
                             "Retentie",
                             "Configuratie"))
Tabel 1.13: Predicties voor random forest
% Voorsp. FALSE % Voorsp. TRUE ID Rij Mtry Min. aantal Retentie Configuratie
53.7% 46.3% validation 968 22 38 TRUE Preprocessor1_Model01
45.4% 54.6% validation 969 22 38 TRUE Preprocessor1_Model01
41.4% 58.6% validation 970 22 38 TRUE Preprocessor1_Model01
41.6% 58.4% validation 971 22 38 TRUE Preprocessor1_Model01
33.4% 66.6% validation 972 22 38 TRUE Preprocessor1_Model01
51.1% 48.9% validation 973 22 38 TRUE Preprocessor1_Model01
21.6% 78.4% validation 974 22 38 FALSE Preprocessor1_Model01
50.0% 50.0% validation 975 22 38 TRUE Preprocessor1_Model01
48.1% 51.9% validation 976 22 38 FALSE Preprocessor1_Model01
22.4% 77.6% validation 977 22 38 TRUE Preprocessor1_Model01
Toon code
## Determine the AUC/ROC curve
rf_auc <- 
  rf_res |> 
  collect_predictions(parameters = rf_best) |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Random Forest")

## Plot the ROC curve
Get_ROC_Plot(rf_auc, position = 2)

## Determine the AUC of the best model
rf_auc_highest   <-
  rf_res |>
  collect_predictions(parameters = rf_best) |> 
  roc_auc(Retentie, .pred_FALSE)

## Add model name and AUC to dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Random Forest", 
          auc = rf_auc_highest$.estimate)
Figuur 1.6: ROC curve voor random forest

1.5 De uiteindelijke fit

  • In de laatste stap van deze analyse maken we het model definitief.
  • We testen het model op de testset en evalueren het model met metrieken en de Variable Importance (VI). De VI kwantificeert de bijdrage van elke variabele aan de voorspellende kracht van een model. Het identificeert welke variabelen significant zijn voor de modelprestaties, wat essentieel is voor het interpreteren en optimaliseren van een model (Van der Laan, 2006). Methoden zoals de Shapley-waarde en permutation importance worden vaak toegepast om dit belang te meten. Op deze methoden komen we terug in het volgende hoofdstuk.

1.5.1 Combineer de AUC/ROC curves en kies het beste model

Eerst combineren we de AUC/ROC curves van de modellen om ze te vergelijken. We kiezen het beste model op basis van de hoogste AUC/ROC.

Toon code
## Combine the AUC/ROC curves to compare the models
Get_ROC_Plot(list(lr_auc, rf_auc))
Figuur 1.7: Gecombineerde ROC curves
Toon code
## Determine which of the models is best based on highest AUC/ROC
dfModel_results <- dfModel_results |>
  mutate(number = row_number()) |> 
  mutate(best = ifelse(auc == max(auc), TRUE, FALSE)) |> 
  arrange(number)

## Determine the best model
sBest_model     <- dfModel_results$model[dfModel_results$best == TRUE]
sBest_model_auc <- round(dfModel_results$auc[dfModel_results$best == TRUE], 4)

Het beste model is het Random Forest model met een AUC/ROC van 0.7. We ronden de analyse verder af met het Random Forest model.

1.5.2 Maak het finale model

We maken het finale model op basis van de beste parameters die we hebben gevonden. Door in de engine bij importance de impurity op te geven, wordt het beste random forest model gekozen om de data definitief mee te classificeren.

## Test the developed model on the test set
## Determine the optimal parameters

## Build the final models
last_lr_mod <-
    logistic_reg(penalty = lr_best$penalty,
                 mixture = 1) |>
    set_engine("glmnet") |>
    set_mode("classification")

last_rf_mod <-
    rand_forest(mtry = rf_best$mtry,
                min_n = rf_best$min_n,
                trees = 1000) |>
    set_engine("ranger", num.threads = cores, importance = "impurity") |>
    set_mode("classification")

1.5.3 Maak de workflow

We voegen het model toe aan de workflow en updaten de workflow met het finale model.

## Update the workflows
 last_lr_workflow <- 
    lr_workflow |> 
    update_model(last_lr_mod)

 last_rf_workflow <- 
    rf_workflow |> 
    update_model(last_rf_mod)

1.5.4 Fit het finale model

We voeren de finale fit uit. De functie last_fit past het model toe op de validatieset.

## Perform the final fit
set.seed(2904)

## Make a final fit for both models so we can save it for later use
last_fit_lr <- 
    last_lr_workflow |> 
    last_fit(splits)

last_fit_rf <- 
    last_rf_workflow |> 
    last_fit(splits)

lLast_fits <- list(last_fit_lr, last_fit_rf) |> 
  set_names(c("Logistic Regression", "Random Forest"))

## Determine which model is best
if(sBest_model == "Logistic Regression") {
  last_fit <- last_fit_lr
} else if(sBest_model == "Random Forest") {
  last_fit <- last_fit_rf
}

## Keep results, model results and associated data
sFittedmodels_outputpath <- Get_Model_Outputpath(mode = "last-fits")
saveRDS(lLast_fits, file = sFittedmodels_outputpath)

sModelresults_outputpath <- Get_Model_Outputpath(mode = "modelresults")
saveRDS(dfModel_results, file = sModelresults_outputpath)

sData_outputpath <- Get_Model_Outputpath(mode = "data")
saveRDS(dfOpleiding_inschrijvingen, file = sData_outputpath)

1.5.5 Evalueer het finale model: metrieken en variable importance

We evalueren het finale model op basis van 4 metrieken: 1) accuraatheid, 2) ROC/AUC en 3) de Brier score (de Mean Squared Error) en 4) de Variable Importance (VI). Uit de VI is op te maken welke variabelen het meest bijdragen aan de voorspelling van de uitkomstvariabele.

## Collect the metrics
last_fit |> 
  collect_metrics() |> 
  mutate(.estimate = round(.estimate, 4)) |>
  knitr::kable(col.names = c("Metriek", 
                             "Estimator",
                             "Estimate",
                             "Configuratie"))
Metriek Estimator Estimate Configuratie
accuracy binary 0.6420 Preprocessor1_Model1
roc_auc binary 0.6491 Preprocessor1_Model1
brier_class binary 0.2219 Preprocessor1_Model1
Toon code
# Extract the variable importance
dfVi <- last_fit |>
  extract_fit_parsnip() |>
  vip::vi() |> 
  arrange(desc(Importance)) |>
  head(20)
  
# Create the plot with fill on the variable 'Importance'
importance_plot <- dfVi |> 
  ggplot(aes(x = reorder(Variable, Importance), 
             y = Importance, 
             fill = Importance)) +
  geom_col(show.legend = FALSE) +
  
  ## Create the title and caption
  labs(x = NULL,
       y = "VI-score",
       caption = sCaption) +
  
  theme_minimal() +
  Set_LTA_Theme() +
  
  theme(
    axis.title.x = element_text(margin = margin(t = 20))
  ) +
  
  coord_flip()
  
  ## Add LTA elements
  importance_plot <- Add_LTA_Theme_Elements(importance_plot, title_subtitle = TRUE)

# Show the plot
print(importance_plot)
Figuur 1.8: Meest voorspellende factoren op basis van de Variable Importance (VI)

1.5.6 Plot de ROC curve

Tot slot visualiseren we de prestaties weer met een ROC curve van het beste model.

## Show the roc curve
auc_lf <- last_fit |> 
  collect_predictions() |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Last fit")

Get_ROC_Plot(auc_lf, position = 3)
Figuur 1.9: ROC curve finale model

1.6 Conclusies

1.6.1 Het beste prognosemodel voor deze opleiding

Het beste prognosemodel blijkt het Random Forest model te zijn.

  • Van de prognosemodellen die we hebben ontwikkeld om retentie na 1 jaar te voorspellen, had het Random Forest model de hoogste AUC/ROC waarde (0.7).

1.6.2 Mate van accuraatheid en lift

Een prognosemodel moet minimaal beter presteren dan een base-model om waarde op basis van accuraatheid toe te voegen. Het base-model neemt als basis de grootste klasse van de gemiddelde retentie na 1 jaar van de afgelopen jaren. Stel we zouden tegen alle studenten zeggen dat ze hun studie gaan halen, dan is de mate van accuratesse gelijk aan dit base-model. Dit base-model is dus altijd hoger dan de 50% lijn van de AUC/ROC curve, tenzij het base-model toevallig precies 50% is.

Toon code
knitr::include_graphics(here::here("R/images", "basemodel-lift.png"))
Figuur 1.10: Lift afhankelijk van base-model en accuraatheid

De mate van accuraatheid van het prognosemodel is vrij laag (64.2%).

  • Base-model: 62.12% – Voor deze opleiding berekenen we het base-model als volgt. Van alle studenten studeerde 62.12% door; 100% - 62.12% = 37.88% studeerde niet door. De grootste klasse van deze twee, 62.12%, is daarmee de accuratesse van het base-model.
  • Accuratesse prognose: 64.2% – Het model voorspelt Retentie na 1 jaar met een accuratesse van 64.2%.
  • Lift: 2.08% – Het model scoort in de huidige opbouw met een verschil van 2.08% (de lift) iets beter dan de accuraatheid van het base-model.

1.6.3 Confusion Matrix

Toon code
## Determine the confusion matrix
confusion_matrix <- last_fit |>
  collect_predictions() |>
  conf_mat(truth = Retentie, estimate = .pred_class) 

dfConf_matrix <- as_tibble(confusion_matrix$table) |>
  rename(Werkelijkheid = Truth) |>
  mutate(Werkelijkheid = ifelse(Werkelijkheid == "TRUE", "Retentie", "Geen retentie"),
         Prediction    = ifelse(Prediction == "TRUE", "Retentie", "Geen retentie"))

pTP  <- Change_Number_Marks((dfConf_matrix$n[4]/sum(dfConf_matrix$n)*100),1)
pFP  <- Change_Number_Marks((dfConf_matrix$n[2]/sum(dfConf_matrix$n)*100),1)
pTN  <- Change_Number_Marks((dfConf_matrix$n[1]/sum(dfConf_matrix$n)*100),1)
pFN  <- Change_Number_Marks((dfConf_matrix$n[3]/sum(dfConf_matrix$n)*100),1)
pACC <- Change_Number_Marks(Last_fit_Accuracy,1)

De prestaties van het model kunnen we verder uitdrukken in een confusion matrix. Hierin zien we de voorspellingen van het model en de werkelijke uitkomsten. De matrix geeft inzicht in de mate van correcte en incorrecte voorspellingen. Ter illustratie werken we de matrix uit voor een voorspelling waarop een bindend studieadvies (BSA) gebaseerd zou kunnen zijn.

Toon code
knitr::include_graphics(here::here("R/images", "confusion-matrix-retention-lta-hhs.png"))
Figuur 1.11: Confusion matrix in relatie tot BSA

We passen de confusion matrix nu toe op het model dat als beste naar voren kwam. De accuraatheid van dit model is 64,2%. De accuraatheid van het model berekenen we door de som van de diagonaal te berekenen: het aandeel goed voorspelde uitkomsten, Retentie = Retentie (True Positive) en Geen retentie = Geen retentie (True Negative), af te zetten tegen het totaal aantal voorspellingen: 57,4% + 6,8% = 64,2%. (NB. De weergave in deze confusion matrix is diagonaal gespiegeld vergeleken met het voorbeeld.)

Toon code
confusion_plot <- plot_confusion_matrix(
    dfConf_matrix,
    target_col = "Werkelijkheid",
    prediction_col = "Prediction",
    counts_col = "n",
    palette = "Blues",
    add_sums = TRUE,
    theme_fn = ggplot2::theme_light,
    sums_settings = sum_tile_settings(
      palette = "Greens",
      label = "Totaal",
      tc_tile_border_color = "black"
    )) +
    
    ## Customize the labels
    labs(
      x = "Werkelijke uitkost",
      y = "Voorspelde uitkomst",
      caption = sCaption
    ) +
    
    Set_LTA_Theme()
  
  ## Add LTA elements
  confusion_plot <- Add_LTA_Theme_Elements(confusion_plot, 
                                           title_subtitle = TRUE)
  
  print(confusion_plot)
Figuur 1.12: Confusion matrix ten opzichte van Retentie na 1 jaar

1.6.4 Uitleggen of verklaren?

Naast de accuraatheid van het model is het ook belangrijk om te weten welke factoren het meest bijdragen aan de voorspelling van retentie na 1 jaar. Daarin gaat de vergelijking met de prestaties van het basemodel mank. Dat model geeft op geen enkele manier aan waarom een student een kans op succes heeft, anders dan - ‘dit is gebruikelijk in deze opleiding’.

Ongeacht de mate van accuraatheid, is het voor onderzoek naar kansengelijkheid essentieel om te weten welke factoren het meest bijdragen aan de voorspelling van retentie na 1 jaar. Het gaat erom dat we het belang van de factoren in de voorspellingen kunnen begrijpen en duiden. Machine Learning is hiervoor uitstekend geschikt, omdat het de mogelijkheid biedt om de belangrijkste factoren en hun invloed te leren kennen (Shmueli, 2010; Shmueli & Koppius, 2011).

1.7 Vervolgstappen: Factoranalyse

De volgende stap (stap 2) is een verdiepende analyse van de mate waarin de factoren die we gevonden hebben van invloed zijn op Retentie na 1 jaar. We kijken naar de rangorde, of ze retentie na 1 jaar verhogen of juist verlagen en hoe stabiel de factoren zijn als we in andere volgordes aan het model toevoegen. Om het concreet te maken zullen we het model toepassen op een aantal fictieve studenten, die we opbouwen uit de meeste voorkomende waarden in deze opleiding. Dit is het onderwerp van analyse 2: de Factoranalyse.

 

Verantwoording

Deze analyse maakt deel uit van het onderzoek naar kansengelijkheid van het lectoraat Learning Technology & Analytics van De Haagse Hogeschool: No Fairness without Awareness | Het rapport is door het lectoraat ontwikkeld in Quarto 1.6.39. | Template versie: 0.9.1.9000

 

Copyright

Dr. Theo Bakker, Lectoraat Learning Technology & Analytics, De Haagse Hogeschool © 2023-2025. Alle rechten voorbehouden.